plot.Abundance<- function (data=NULL,station="ALL", species="ALL", year="ALL", group.species=FALSE, group.station=FALSE,...) {
    abund.data<-data
    if (is.null(abund.data)==TRUE) stop("A data file is required.")
    abund.data$seAbund<-ifelse(is.na(abund.data$seAbund)==TRUE,0,abund.data$seAbund)
    abund.data$lower95<-abund.data$MeanAbund-2*abund.data$seAbund
    abund.data$upper95<-abund.data$MeanAbund+2*abund.data$seAbund
    abund.data$lowerSE<-abund.data$MeanAbund-abund.data$seAbund
    abund.data$upperSE<-abund.data$MeanAbund+abund.data$seAbund
    station.list<-station
    if (station[1]=="ALL") station.list<-levels(factor(abund.data$station_name)) # really need to make this station!!!
    species.list<-species
    if (species[1]=="ALL") species.list<-levels(factor(abund.data$species))
    if (group.species==FALSE & group.station==FALSE) {
        for (Species in species.list) {
            for (Station in station.list) {
                dataset<-abund.data[abund.data$species==Species & abund.data$station_name==Station,]
                dataset$year<-factor(as.character(dataset$year))
                dataset<-dataset[order(dataset$year),]
                if (nrow(dataset)>0 & suppressWarnings(max(dataset$upperSE>0))) {        
                    trellis.device(device=windows,new=FALSE, col=FALSE) 
                    print(xYplot(Cbind(MeanAbund,lowerSE,upperSE)~as.integer(as.character(year)),
                        data=dataset,
                        type='b',
                        scales = list(x=list(rot=45,at=seq(min(as.integer(as.character(dataset$year))),length=max(as.integer(as.character(dataset$year)))-min(as.integer(as.character(dataset$year)))+1)),cex=1.5,cex.axis=2,cex.main=2.25),
                        xlab = list(label="Year",cex=1.5),
                        #xlim=(c(min(as.integer(as.character(dataset$year)))-1,max(as.integer(as.character(dataset$year)))+1)),
                        ylim=(c(max(c(0,min(dataset$lowerSE)-0.1*min(dataset$lowerSE))),1.1*max(dataset$upperSE))),
                        pch=19,
                        cex=1.5,
                        lty=1,
                        lwd=2,
                        label.curves=FALSE,
                        ylab=list(label="Numbers of birds detected",cex=1.5),
                        main=list(label=paste("Detection Estimates \n for",dataset[1,]$species,"\n at",dataset[1,]$station_name),cex=1.6),
                        )
                    )
                }
            }
        }
    }
    if (group.species==TRUE & group.station==FALSE) {
        for (Station in station) {
            dataset<-abund.data[abund.data$species %in% species.list & abund.data$station_name==Station,]
            dataset$species<-factor(dataset$species)
            dataset$year<-factor(as.character(dataset$year))
            dataset<-dataset[order(dataset$year),]
            if (nrow(dataset)>0 & suppressWarnings(max(dataset$upperSE>0))) {        
                trellis.device(new=FALSE, col=FALSE) 
                print(dataset=dataset,x=xYplot(Cbind(MeanAbund,lowerSE,upperSE)~as.integer(as.character(year)),
                    data=dataset,
                    groups=species,
                    type='b',
                    scales = list(x=list(rot=45),cex=1.5,cex.axis=2,cex.main=2.25,tick.number=length(min(as.integer(as.character(dataset$year))):max(as.integer(as.character(dataset$year))))),
                    xlab = list(label="Year",cex=1.5),
                    xlim=(c(min(as.integer(as.character(dataset$year)))-1,max(as.integer(as.character(dataset$year)))+1)),
                    ylim=(c(max(c(0,min(dataset$lowerSE)-0.1*min(dataset$lowerSE))),1.1*max(dataset$upperSE))),
                    cex=1.5,
                    lwd=2,
                    keys='lines',
                    ylab=list(label="Number of birds detected",cex=1.5),
                    main=list(label=paste("Abundance Indices \n at",dataset[1,]$station_name),cex=1.6)
                    )
                )
            }
        }
    }
    if (group.species==FALSE & group.station==TRUE) {
        for (Species in species.list) {
            dataset<-abund.data[abund.data$species==Species & abund.data$station_name %in% station.list,]
            dataset$station_name<-factor(dataset$station_name)
            dataset$year<-factor(as.character(dataset$year))
            dataset<-dataset[order(dataset$year),]
            if (nrow(dataset)>0 & suppressWarnings(max(dataset$upperSE>0))) {        
                trellis.device(new=FALSE, col=FALSE) 
                print(dataset=dataset,x=xYplot(Cbind(MeanAbund,lowerSE,upperSE)~as.integer(as.character(year)),
                    data=dataset,
                    groups=station_name,
                    type='b',
                    scales = list(x=list(rot=45),cex=1.5,cex.axis=2,cex.main=2.25,tick.number=length(min(as.integer(as.character(dataset$year))):max(as.integer(as.character(dataset$year))))),
                    xlab = list(label="Year",cex=1.5),
                    xlim=(c(min(as.integer(as.character(dataset$year)))-1,max(as.integer(as.character(dataset$year)))+1)),
                    ylim=(c(max(c(0,min(dataset$lowerSE)-0.1*min(dataset$lowerSE))),1.1*max(dataset$upperSE))),
                    cex=1.5,
                    lwd=2,
                    keys='lines',
                    ylab=list(label="Numbers of birds detected",cex=1.5),
                    main=list(label=paste("Abundance Indices \n for",dataset[1,]$species),cex=1.6)
                    )
                )
            }
        }
    }    
}
